Example Usage

[2]:

import chalc.chromatic as chromatic
from chalc.sixpack import SubChromaticInclusion, SixPack
from chalc.plotting import plot_sixpack, plot_diagram, draw_filtration, animate_filtration
import numpy as np
import matplotlib.pyplot as plt

For our data we sample 100 points on a circle with some noise and 100 points from inside the unit disk.

[4]:

rng = np.random.default_rng(40)
num_points_circle = 100
num_points_disk = 100
mean = [0, 0]
cov = np.eye(2)*0.01
circle = np.array([[np.sin(2*np.pi*t), np.cos(2*np.pi*t)] for t in rng.random(num_points_circle)]).T +\
    rng.multivariate_normal(mean, cov, num_points_circle).T # points as columns
disk = np.sqrt(rng.random(num_points_disk)) * np.array([[np.sin(2*np.pi*t), np.cos(2*np.pi*t)] for t in rng.random(num_points_disk)]).T
points = np.concatenate((circle, disk), axis=1)
plt.scatter(circle[0, :], circle[1, :], s=10)
plt.scatter(disk[0, :], disk[1, :], s=10)
plt.gca().set_aspect('equal')
colours = [0]*num_points_circle + [1]*num_points_disk
plt.title('Point cloud')
plt.show()
_images/example_5_0.svg

We can compute the chromatic alpha/Delaunay–Čech/Delaunay–Rips complex filt of the point cloud using the functions chromatic.alpha, chromatic.delcech, and chromatic.delrips respectively. In each of these cases, filt has far fewer simplices than either the Čech or Vietoris–Rips complex, which each have \(\displaystyle \binom{200}{2} = 19900\) edges and \(\displaystyle \binom{200}{3} = 1313400\) 2-simplices.

[5]:

# max_num_threads=0 gives automatic parallelism
chromatic_alpha_filt, numerical_instability_flag = chromatic.alpha(points, colours, max_num_threads=0)
print(f'\'chromatic_alpha_filt\' has {len(chromatic_alpha_filt.simplices[1])} 1-simplices and {len(chromatic_alpha_filt.simplices[2])} 2-simplices.')

'chromatic_alpha_filt' has 954 1-simplices and 1312 2-simplices.

chromatic_alpha_filt is an object of type Filtration, and simplices in the filtration are objects of type Simplex. Both of these have convenience methods:

[6]:

# __len__()
print(f'\'chromatic_alpha_filt\' has {len(chromatic_alpha_filt)} simplices and {chromatic_alpha_filt.num_vertices} vertices.')
print(f"Dimension of the filtration: {chromatic_alpha_filt.dimension}.")

facets = []
# __iter__()
for simplex in chromatic_alpha_filt:
    # Do something with the simplex
    if simplex.dimension == 2:
        print(f"Found a simplex of dimension {simplex.dimension} with vertices {simplex.vertices} and colours {simplex.colours}")
        # Get a list of handles to the facets of the simplex
        facets = simplex.facets
        break

if facets:
    verts = facets[0].vertices
    # __contains__()
    if verts in chromatic_alpha_filt:
        print(f"{verts} is a simplex in the filtration.")
    else:
        print(f"{verts} is not a simplex in the filtration.")

'chromatic_alpha_filt' has 3023 simplices and 200 vertices.
Dimension of the filtration: 3.
Found a simplex of dimension 2 with vertices [30, 33, 87] and colours [0]
[33, 87] is a simplex in the filtration.

We can visualize the filtration at any given time using plotting.draw_filtration. For example, at time 0.3 we have:

[7]:

fig, ax = plt.subplots(1, 3, figsize=(12, 4))
draw_filtration(chromatic_alpha_filt, points, time=0.3, include_colours=[0], ax=ax[0])
ax[0].set_title('Colour 0')
draw_filtration(chromatic_alpha_filt, points, time=0.3, include_colours=[1], ax=ax[1])
ax[1].set_title('Colour 1')
draw_filtration(chromatic_alpha_filt, points, time=0.3, ax=ax[2])
ax[2].set_title('Both colours')
for i in range(3):
    ax[i].set_xlim(-1.2, 1.2)
    ax[i].set_ylim(-1.2, 1.2)
fig.suptitle('Complexes at $t=0.3$')
plt.tight_layout()
plt.show()
_images/example_11_0.svg

Next we compute the 6-pack of persistent homology diagrams of the chromatic alpha/Delaunay–Čech complex/Delaunay–Rips complexes. To do this we need to choose a map of filtrations induced by some inclusion of colours. For this example, we will consider the inclusion of the blue points (colour 0) into all points by using the SubChromaticInclusion class.

[8]:

# Using the previously computed filtration:
dgms_alpha = SubChromaticInclusion(chromatic_alpha_filt, {0}).sixpack()
# Using the Delaunay--Cech and Delaunay--Rips filtrations
chromatic_delcech_filt, _ = chromatic.delcech(points, colours)
chromatic_delrips_filt, _ = chromatic.delrips(points, colours)
dgms_delcech = SubChromaticInclusion(chromatic_delcech_filt, {0}).sixpack()
dgms_delrips = SubChromaticInclusion(chromatic_delrips_filt, {0}).sixpack()

In code above, SubChromaticInclusion corresponds to the inclusion of some subset (or subsets) of colours into the whole filtration. If we had more colours, say {0, 1, 2, 3}, we might also consider the inclusion of the sub-filtration spanned by simplices whose colours lie in any subset among (say) {{0, 1}, {1, 3}, {2}}. This would be specified by SubChromaticInclusion([[0, 1], [1, 3], [2]]) above. You can also consider other mapping methods; see KChromaticInclusion, SubChromaticQuotient and KChromaticQuotient.

Each 6-pack is an object of type SixPack that has the kernel, cokernel, domain, codomain, and relative persistent homology diagrams. The class SixPack implements the abstract base class Mapping from the Python standard library and supports dictionary-like syntax. In particular, you can access the diagrams by indexing into the 6-pack.

[9]:

print(f"The diagram names are: {list(dgms_alpha.keys())}")

# Iterate over the diagrams like in a dictionary.
for diagram_name in dgms_alpha:
    print(diagram_name)
    # do something with dgms_alpha[diagram_name]

# Access individual diagrams by indexing.
ker = dgms_alpha['ker']

The diagram names are: ['ker', 'cok', 'dom', 'cod', 'im', 'rel']
ker
cok
dom
cod
im
rel

Every feature in the diagram is represented by either a pair of simplices \((\sigma, \tau)\) if the feature has a finite birth and death time, or an unpaired simplex \(\sigma\) if the feature does not die. For a \(p\)-dimensional feature, \(\dim(\sigma) = p+1\) if the diagram is the kernel, and \(\dim(\sigma) = p\) otherwise.

[10]:

print(dgms_alpha["ker"])

Paired: frozenset({(379, 737), (367, 418), (289, 358), (996, 1010), (494, 623), (228, 235), (2118, 2957), (1498, 1504), (562, 568), (390, 477), (1080, 1110), (507, 644), (479, 506), (618, 707), (1035, 1047), (828, 919), (392, 553), (1122, 1143), (742, 986), (550, 584), (349, 484), (1403, 1406), (1295, 1446), (653, 862), (1133, 1139), (519, 1023), (254, 521), (299, 428), (2084, 2091), (264, 412), (884, 938), (401, 433), (814, 1020), (832, 1754), (372, 702), (743, 746), (438, 842), (886, 907), (1245, 1310), (287, 328), (1013, 1018), (1489, 1490), (1029, 1096), (351, 825), (999, 1098), (608, 609), (389, 470)})
Unpaired: frozenset()

Each individual diagram is an object of type SimplexPairings, which supports set-like syntax and implements the abstract base class Collection from the Python standard library.

[11]:

# __bool__()
print("The kernel is non-empty." if dgms_alpha["ker"] else "The kernel is empty.")
# Or you can use bool(dgms_alpha.ker).

# __len__()
# Get the number of features in the diagram
print(f"The cokernel has {len(dgms_alpha['cok'])} features.")

# __contains__()
# Check if a feature exists in a diagram.
print("The domain contains a feature represented by the simplex pair (20, 40): "
      f"{(20, 40) in dgms_alpha["dom"]}")

# __iter__()
# Iterate over the paired and unpaired simplices.
for feature in dgms_alpha['ker']:
    if isinstance(feature, tuple):
        sigma, tau = feature
        # Do something with a pair of simplices.
        ...
    else:
        sigma = feature
        # Do something with an unpaired simplex.
        ...

# The paired and unpaired simplices can be considered separately.
print("The \"paired\" property of the kernel has type: "
      f"{type(dgms_alpha['ker'].paired)}")

print("The \"unpaired\" property of the kernel has type: "
      f"{type(dgms_alpha['ker'].unpaired)}")

for sigma, tau in dgms_alpha['ker'].paired:
    # Do something with paired simplices.
    ...

for sigma in dgms_alpha['dom'].unpaired:
    # Do something with unpaired simplices.
    ...

The kernel is non-empty.
The cokernel has 234 features.
The domain contains a feature represented by the simplex pair (20, 40): False
The "paired" property of the kernel has type: <class 'frozenset'>
The "unpaired" property of the kernel has type: <class 'frozenset'>

A convenience function get_matrix is provided to extract the diagrams as a matrix of birth/death times:

[12]:

np.set_printoptions(threshold=10)
# Get the domain in dimension zero.
print(dgms_alpha.get_matrix('dom', 0))
# Get the kernel in dimensions zero and one.
print(dgms_alpha.get('ker', [0, 1]))
# Get the kernel in all dimensions from zero to max(dgms_alpha.dimensions).
print(dgms_alpha.get('ker'))

[[0.         0.01367401]
 [0.         0.01173629]
 [0.         0.02611904]
 ...
 [0.         0.04149252]
 [0.         0.06114399]
 [0.                inf]]
Paired: frozenset({(379, 737), (367, 418), (289, 358), (996, 1010), (494, 623), (228, 235), (2118, 2957), (1498, 1504), (562, 568), (390, 477), (1080, 1110), (507, 644), (479, 506), (618, 707), (1035, 1047), (828, 919), (392, 553), (1122, 1143), (742, 986), (550, 584), (349, 484), (1403, 1406), (1295, 1446), (653, 862), (1133, 1139), (519, 1023), (254, 521), (299, 428), (2084, 2091), (264, 412), (884, 938), (401, 433), (814, 1020), (832, 1754), (372, 702), (743, 746), (438, 842), (886, 907), (1245, 1310), (287, 328), (1013, 1018), (1489, 1490), (1029, 1096), (351, 825), (999, 1098), (608, 609), (389, 470)})
Unpaired: frozenset()
Paired: frozenset({(379, 737), (367, 418), (289, 358), (996, 1010), (494, 623), (228, 235), (2118, 2957), (1498, 1504), (562, 568), (390, 477), (1080, 1110), (507, 644), (479, 506), (618, 707), (1035, 1047), (828, 919), (392, 553), (1122, 1143), (742, 986), (550, 584), (349, 484), (1403, 1406), (1295, 1446), (653, 862), (1133, 1139), (519, 1023), (254, 521), (299, 428), (2084, 2091), (264, 412), (884, 938), (401, 433), (814, 1020), (832, 1754), (372, 702), (743, 746), (438, 842), (886, 907), (1245, 1310), (287, 328), (1013, 1018), (1489, 1490), (1029, 1096), (351, 825), (999, 1098), (608, 609), (389, 470)})
Unpaired: frozenset()

You can also access the entrance_times and dimensions of simplices in the 6-pack diagrams.

[13]:

print(dgms_alpha.entrance_times)
print(dgms_alpha.dimensions)

[ 0.          0.          0.         ... 16.97606576 16.97606576
 16.97606576]
[0 0 0 ... 2 2 3]

The convenience functions plotting.plot_sixpack and plotting.plot_diagram are provided for plotting either a 6-pack or a specific diagram from a 6-pack.

[14]:

fig1, ax1 = plot_sixpack(dgms_alpha, threshold = 1e-3)
fig1.suptitle('Chromatic Alpha')
fig1.set_figwidth(15)
fig1.set_figheight(9)
fig1.subplots_adjust(top=0.92)
plt.show()
fig2, ax2 = plt.subplots(1, 3)
# You can specify the dimensions of features to include in the plot.
plot_diagram(dgms_alpha, 'rel', dimensions = {0, 2}, truncation = 0.3, ax = ax2[0], threshold = 1e-3)
ax2[0].set_title('Chromatic Alpha')
# You can also specify a single dimension as an integer.
plot_diagram(dgms_delcech, 'rel', dimensions = 1, truncation = 0.3, ax = ax2[1], threshold = 1e-3)
ax2[1].set_title('Chromatic Delaunay-Cech')
# If no dimensions are specified, all features will be included.
plot_diagram(dgms_delrips, 'rel', truncation = 0.3, ax = ax2[2], threshold = 1e-3)
ax2[2].set_title('Chromatic Delaunay-Rips')
fig2.suptitle('Individual relative diagrams')
fig2.set_figwidth(15)
plt.show()
_images/example_25_0.svg
_images/example_25_1.svg

We can save the diagrams to a HDF5 file or load a diagram from a HDF5 file. The HDF5 format is

[15]:

with h5py.File('test.h5', 'w') as f:
    dgms_alpha.save(f)

with h5py.File('test.h5', 'r') as f:
    dgms_alpha_from_file = SixPack.from_file(f)

print(dgms_alpha_from_file == dgms_alpha)

True

We can visualize the 2-skeleton of a chromatic filtration for points in 2D:

[16]:

animation = animate_filtration(
    chromatic_alpha_filt, points, filtration_times=np.linspace(0, 1.0, 45).tolist(), animation_length=5)
animation

[16]:
_images/example_29_1.svg